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ABSTRACT 

X 

This paper presents a fully three-dimensional radiative hydrodymanics simu- 
lation with realistic opacities for a gravitationally unstable 0.07 M disk around 
a 0.5 M star. We address the following aspects of disk evolution: the strength of 
gravitational instabilities under realistic cooling, mass transport in the disk that 
arises from GIs, comparisons between the gravitational and Reynolds stresses 
measured in the disk and those expected in an a-disk, and comparisons between 
the SED derived for the disk and SEDs derived from observationally determined 
parameters. The mass transport in this disk is dominated by global modes, and 
the cooling times are too long to permit fragmentation for all radii. Moreover, 
our results suggest a plausible explanation for the FU Ori outburst phenomenon. 



Subject headings: accretion, accretion disks — convection — hydrodynamics - 
instabilities — solar system: formation 
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1. 



INTRODUCTION 



Accurately simulating cooling in protoplanetary disks is fundamental to modeling disk 
evolution. If a disk with enough mass becomes sufficiently cold, gravitational instabilities 
(GIs) set in and produce nonaxisymmetric structure, which provides mass and angular mo- 
mentum transport through long-range torques. The strength of these GIs and their efficiency 
in transporting mass depend on the cooling rate (Lodato & Rice 2004; Mejfa et al. 2005). 
Moreover, if the cooling rates are very high, the stress in the disk can become large enough 
to cause fragmentation (Gammie 2001; Rice et al. 2005) and possibly induce formation of 
bound protoplanetary clumps (Boss 1997, 1998, 2001, 2002). Because cooling rates can be 
a principal determinant of disk evolution (Pickett et al. 1998, 2000), we here investigate the 
evolution of a 0.07 M disk surrounding a 0.5 M Q star by using three-dimensional radiative 
hydrodynamics with realistic opacities. 

A gas disk is both massive and cold enough for GIs to set in when the Toomre (1964) 
parameter 



becomes less than about 1.5 (Durisen et al. 2003); here, c s is the sound speed, K e is the 
epicyclic frequency, and £ is the gas surface density. The strength of these GIs strongly 
depends on the thermal physics of the gas (see Pickett et al. 2000) and on disk cooling. To 
date, simulations of gravitationally unstable disks exhibit a variety of evolutionary behav- 
iors. Many of the differences can be linked to the equation of state and the techniques used 
to model radiative cooling (for a review, see Durisen et al. 2006). For example, a global 
cooling time t coo i = constant everywhere is used in the simulations of Pickett et al. (1998, 
2000, 2003, hereafter Paper I) and Mejfa et al. (2005, hereafter Paper II). Their simulations 
show that mass and angular momentum transport are dominated by low-order modes. When 
their disks are in the asympototic phase (see Paper II), where shock heating from gravita- 
tional instabilities balances cooling, a simple a-disk description (Shakura & Sunyaev 1973) 
inaccurately describes the disk evolution. However, if a cooling time that is dependent on 
the the gas orbital speed is enforced, i.e., t coo i^ = constant (e.g., Gammie 2001; Lodato & 
Rice 2004; Mayer et al. 2004), mass and angular momentum transport can be described by 
an a-disk model with an 



as discussed in Gammie (2001). Here, 7' is the two-dimensional ratio of specific heats 
and equals 3 — 2/7 in the strongly self-gravitating limit and (37 — 1) / (7 + 1) in the non- 
self-gravitating limit, where 7 is the three-dimensional ratio of specific heats. A t coo i^ = 




(2) 
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constant leads to an a-disk-like behavior in the asymptotic phase because the stress in the 
disk that results from GIs must lead to heating that balances t coo i^ = constant cooling. 
Similarly, prescribing a global cooling time should lead to a dominance in disk heating 
by lower-order modes, which provide significant long-range torques over most of the disk. 
Let P ro t = rotation period = 2w/Q. For either prescription, when t coo i becomes less than 
about 3/Q ~ -Prot/2, for a 7 = 5/3 gas, the disk fragments (Gammie 2001; Rice et al. 
2003; Mejfa et al. 2005); simulations in Rice et al. (2005) suggest that this critical cooling 
time for fragmentation may be closer to Q/Q rs P rot for 7 = 5/3 and that it increases to 
about 12/Q ~ 2P rot for 7 = 7/5. Because different cooling prescriptions lead to different 
GI behavior and disk evolution, radiative hydrodynamics with realistic opacities, a realistic 
equation of state (EOS), and proper boundary conditions must be used to address how real 
disks will behave. 

Some work has already been done on the protoplanetary disk problem with radiative hy- 
drodynamics. Nelson et al. (2000) used a radiative cooling scheme for their two-dimensional 
hydrodynamics simulations by extrapolating vertical structure from their thin disks by as- 
suming polytropic vertical hydrostatic equilibrium. They also created SEDs from their data 
and compared their disks with SEDs derived from observed parameters. Boss (2001) im- 
plemented fully 3-D radiative diffusion (Bodenheimer et al. 1990). Mejia (2004) employed 
fully 3-D flux-limited diffusion with different boundary conditions. The Boss (2001) and the 
Mejia (2004) results are quite different. Boss (2001, 2002, 2004, 2005) argues that radiative 
physics permits disk fragmentation because convection makes the cooling times very short. 
This fragmentation leads to direct planet formation by disk instability. Furthermore, Boss 
(2002) argues that the disk cooling rate is insensitive to metallicity, and that the observed 
metallicity-planet correlation (Fischer & Valenti 2005), which is often considered support 
for core-accretion plus gas capture (Pollack et al. 1996) when taken at face value, can be 
explained by metallicity's effect on inward planet migration rates (Boss 2005). On the other 
hand, the simulation presented here along with the results of Cai et al. (2006), which both 
use the radiative scheme developed by Mejia (2004), show that cooling times are too long 
to permit direct planet formation by disk fragmentation. However, GIs may still aid planet 
formation. Durisen et al. (2005), following work by Haghighipour & Boss (2003), suggest 
that rings formed by GIs could lead to accelerated core-accretion plus gas capture. In fact, 
Rice et al. (2004) find that solids are concentrated into the spiral arms of GI active disks. 

The subject of Paper II was a series of simulations with a constant global cooling time 
(CCT) prescription. This paper presents a disk simulation evolved under realistic radiative 
cooling with the Mejia (2004) scheme, and it is an extension with new analyses of the work 
presented in Mejia (2004). Section 2 describes the radiative cooling (RC) algorithm used 
in the hydrodynamics code, while §3 presents the simulation setup, evolution, and results. 
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Comparisons between the radiative cooling and constant cooling (Paper II) simulations are 
made in §4. Section 5 discusses the relevance of our work to real disks and compares our 
results to work by other groups. We summarize our conclusions in §6. 
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Fig. 1. — Scheme showing the two optical regions of the initial disk (to scale). The left 
inset represents part of a column where both the atmosphere and the interior of the disk 
are shown. The rectangle with the diagonal pattern is the current atmospheric cell at which 
the cooling is being calculated using equation (3). Optical depths, temperatures, and the 
boundary fluxes are represented by arrows. The right inset shows the boundary conditions 
for the last interior column in the radial direction, where the radial boundary fluxes are set 
equal to the vertical boundary flux. 



2. HYDRODYNAMICS 



The three-dimensional hydrodynamics code with self-gravity is the same as that used 
in Papers I and II and is described in detail in Pickett (1995), Pickett et al. (1998, 2000), 
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and Mejia (2004). The hydrodynamics code solves Poisson's equation, an ideal gas equa- 
tion of state, and the equations of hydrodynamics (Yang 1992) in conservative form on a 
uniform cylindrical grid (r, <p, z). The code computes the source and flux terms (Norman 
& Winkler 1986) separately in an explicit, second-order time integration (van Albada et al. 
1982; Christodoulou 1989; Yang 1992), where the advective terms are calculated with a van 
Leer scheme (Norman & Winkler 1986). The energy equation has the form (Williams 1988; 
Pickett et al. 2000; Mejia 2004) 

a Ah i 

_____ + V ■ (e^v) = -e xh ~ l (r - A - V • F) , (3) 

where e is the internal energy density and the heating term T includes the effects of shock 
heating by artificial bulk viscosity through a second-order Neumann & Richtmeyer scheme 
(see Pickett 1995). This artificial viscosity ensures that the jump conditions are satisfied 
by adding the correct amount of entropy to the gas. For more details on the implemented 
AV scheme, we refer the reader to Pickett (1995) and to Pickett et al. (2000). The cooling 
terms A and V • F, which are described below, represent the local radiative cooling in the 
atmosphere of the disk and in the interior of the disk, respectively. 



2.1. Radiative Cooling 

Our radiative physics scheme employs two approximations for the divergence of the flux 

V-F= / pn(S-J) dtt, (4) 



Itt 



where k is the mass absorption coefficient, p is the density, S is the source function, and 
J is the mean intensity. The two limits are fit together with an Eddington-like boundary 
condition. When the optical depth r is very large, equation (4) takes a diffusion form, which 
is approximated in finite difference form as 

1 A (rF r ) 1 AFs AF Z 

V-F= )-Jl + -—A + —JL. 5 

r Ar r A<p Az 

In this limit, the radiative flux is calculated by flux-limited diffusion, as described below, 
and Rosseland mean opacities, where k = XRoss (see Appendix A), are used. When the mean 
intensity becomes negligible and a parcel of gas is allowed to radiate as much as its emissivity 
allows, equation (4) takes the form 



V • F = 4pK Planck aT 4 , 



(6) 
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where the Planck mean opacity ftpianck is appropriate in this limit and the emissivity B(T)Kpi anc k 

(aT 4 /ll) Kpi anck . 

We define two regions of the disk, the interior and the atmosphere, through the Rosse- 
land optical depth, which is integrated vertically downward toward the disk midplane. The 
interior is defined for cells with tr oss > 2/3, and in this region the cooling algorithm uses 
flux-limited diffusion. When rp oss < 2/3, which defines the atmosphere cells, a different ap- 
proach is required. A real atmosphere would not cool according to equation (6), but it would 
absorb photons emerging from the interior of the disk. Therefore, we couple the atmosphere 
to the interior by defining the net atmosphere cooling as 

A = 4pK Planck aT 4 - pXRossF b dr y e" rup , (7) 

where the first term on the right-hand side indicates pure radiative cooling for an atmospheric 
cell given by its own temperature and the second term on the right-hand side accounts for 
heat gained by the atmospheric cell from the underlying disk. The second term is nonzero if 
there is an upward energy flux Fbdry at the boundary between the atmosphere and the interior 
of the disk. The boundary flux is diminished exponentially by the Rosseland optical depth, 
which is measured from the upper boundary cell (labeled r up in Fig. 1). If the boundary 
flux is negative (downward) due to shocks and/or irradiation in the atmosphere, the second 
term on the right in equation (7) is made zero such that the cell can cool as much as its 
emissivity allows. Without the second term on the right-hand side of (7), our atmosphere 
regions contract into a single cell within a few outer rotations. 

The boundary flux provides the boundary condition between the atmosphere and the 
interior of the disk by defining the balance between the energy leaving the disk interior 
through the atmosphere and the energy gained from the atmosphere. If no atmospheric 
heating is assumed and all the energy that flows through the boundary is radiated with 
an Eddington effective temperature fitted to some boundary temperature Tbd ry and the 
corresponding Rosseland optical depth at Tbdry, both of which are measured at the center of 
the first disk interior cell in a column, the flux at the boundary would be 



bdry 



3 (r bdry + 2/3) 



similar to the equation for Eddington gray atmosphere. To add the contribution from the 
atmosphere, it is assumed that half the radiation emitted by the atmosphere leaves the top 
of the atmosphere and half shines on the disk. Then the boundary flux becomes 

p _ 4(7 ( T bdry ~ T atm) (Q , 

bdry " 3(r bdry + 2/3) ' W 
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where the atmospheric T atm is a measure of the downward flux into the interior of the disk 
from the atmosphere, and it is given by half the sum of the atmospheric cooling 

^a 4 t m = 2^pK Planck T 4 A^ (10) 

atm 

Equation (10) only accounts for the first term on the right in equation (7). In the 3-D code, 
^atm is calculated first, then Fbd ry , and finally A. Because the atmospheric cells can be heated 
by shocks generated by Gl-activity (Boley & Durisen 2006) and because these routines are 
also used to treat outside irradiation in other calculations (Cai et al. 2006), it is possible for 
^atm > ^bdry which makes Fbdry < 0. These somewhat cumbersome procedures are necessary 
to fit the explicitly modeled upper, 1ow-tr oss region to the inner, high-ra oss region through 
an Eddington-like solution across the bounding cell, where the Eddington-like solution fitted 
across this cell has to be modified, as in equation (9), to include downward radiation from 
above the cell. 

The discussion above only applies to the vertical boundary condition in grid columns 
that contain both atmosphere and interior cells, like the first inset in Figure 1. However, 
another boundary condition must be set to the fluxes between cells when one of them is 
in the interior and an adjacent one, in r or <p, is in the atmosphere. Radial and azimuthal 
boundary fluxes are uncoupled with the adjacent atmospheric cells because, in effect, A 
is considered to act only in the vertical direction. Therefore, the values of the r and 
boundary fluxes are made equal to the vertical boundary fluxes under the assumption that 
the atmosphere surrounding an interior column has about the same properties in all other 
directions as it does in z. Hence, a boundary cell will radiate (or gain) energy through a face 
contiguous with an atmosphere cell in the same manner as the column's vertical boundary. 
This crude stair-step model for the disk surface enlarges the surface area and probably errs 
in the direction of enhanced disk cooling rates. 

Finally, for the interior of the disk, the flow of energy is calculated by flux-limited 
diffusion, i.e., 

XrossP Ax 

for radiative flux F x in direction x = (r, <fr, z) and Ax = (Ar, A<fr, Az) on a cylindrical grid. 
Because fluxes are measured at cell faces, XRoss, P, and T in equation (11) are the averages of 
the two adjacent cells that share the face at which the flux is being calculated. The quantity 
(3 is the flux limiter that, according to Bodenheimer et al. (1990), has the form 



= 2 + y 

6 + 3y + y 2 ' 



(12) 
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for 

4 

y = 7f 

The flux limiter is 1/3 in the high optical depth limit, while, in the optically thin limit, (3 
asymptotes to 1/y, making F x — > 4crT 4 . Remember, however, that the diffusion approxima- 
tion is only used up to the boundary cell where tr oss > 2/3. 

Equation (11) is used to calculate all fluxes in the interior of the disk with four excep- 
tions. (1) When an adjacent cell is in the atmosphere, the flux at the common face between 
both cells is set to F bdry as discussed above. (2) If the atmosphere extends to the midplane, 
Fbdry = 0. (3) In the case that the interior extends all the way to the top of the grid, the flux 
at the top face is given by equation (8). (4) All vertical fluxes at the equatorial plane are 
because of the mirror symmetry assumption. In the disk interior, we use the flux determined 
by equation (11) in equation (5) to calculate the divergence of the flux; for this reason, we 
refer to the cooling in the disk interior simply as V • F cooling. Appendix B demonstrates 
the accuracy of this algorithm. 

Because the local radiative time scale, as defined by 

e e . . 

t A or t V F = ^ or — — , (14a) 

can be much shorter than the Courant condition, we limit the local cooling time to be no less 
than about 10% of the initial outer rotation period (ORP) of the gas disk at about 33 AU, 
where 1 ORP ~ 250 yr. Without this limiting, saw-toothing in the temperature structure 
of the disk becomes noticeable, and the At becomes small and computationally prohibitive. 
For similar reasons, the local heating time scale by artificial viscosity 

*AV = TT (14b) 
J- AV 

is also limited to be no less than about 0.1 ORPs. Preliminary tests without these limits 
show that cooling and viscous heating times rarely become less than 0.1 ORPs when the 
disk has evolved away from its initial state and the cooling/heating-limited cells tend to 
be confined to uninteresting parts of the calculation. These local time scale limiters will 
not prevent fragmentation, because fragmentation occurs for i cool « P TOt for a 7 = 5/3 disk 
(Gammie 2001; Rice et al. 2003; Mejfa et al. 2005; Rice et al. 2005). 



AT 



Ax 



(13) 
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3. THE SIMULATION 

3.1. Initial Conditions 

The initial model is the same one used for the CCT simulations described in §3.1 of 
Paper II, i.e., a nearly Keplerian, equilibrium disk of 0.07 M that extends from 2.3 to 
40.0 AU and orbits a star of 0.5 M Q . The initial surface density profile is S ~ r -°- 5 an d the 
initial equatorial temperature profile is T ~ r" 1 over most of the radial range. Opacities from 
D'Alessio et al. (2001) are used to compute Xross and Kpianck (see Appendix A). Temperatures 
are calculated every time step by assuming an ideal gas law, with densities and pressures 
already calculated in the code, and by iteratively solving for the molecular weight (table taken 
from D'Alessio 1996), because it too is a function of pressure or density and temperature. Due 
to an error with the inclusion of helium in the solar mix, the mean molecular weight used for 
the temperature and pressure ranges in the simulation is around 2.7 when it should be close 
to 2.3. This introduces a systematic offset no larger than about 16% into the temperature, 
which directly affects the cooling rates. However, as discussed below, the opacity law is 
roughly quadratic in temperature, which means the flux calculated from equation (11) is 
roughly quadratic in temperature, too. The error in the mean molecular weight should be 
an error that typically enhances the cooling by 1.16 2 or about 4/3. Because we expect the 
cooling to be enhanced, fragmentation should be more likely in the RC disk than it would 
be with the correct mean molecular weights. 

The minimum allowed temperature at all times is 3 K to simulate radiating into empty 
space. The temperature of this disk never reaches the dust sublimation temperature, near 
T ~ 1400 K (Muzerolle et al. 2003), so the opacity is mostly due to dust. A maximum grain 
size a max = 1/xm (see Appendix A) is used. 

The main simulation presented in this paper evolves for a little over 16 ORPs, or about 
4,000 yr. The resolution used is the same as all the constant cooling time simulations of 
Paper II, but the z direction is fixed at 32 zones. The initial disk extends almost to the edge 
of the starting grid, namely (r,<p,z) = (256,128,32). The r, z cross-section of this disk in 
the initial grid is accurately portrayed in Figure 1. The grid is extended to 512 in r when 
the disk expands after the first few ORPs. Two phases of the simulation are also run in a 
high azimuthal resolution grid, (r, <fr, z) = (512,512,32) to test for fragmentation (see §5.2). 
A random cell-to-cell density perturbation of amplitude | Ap/ p\ = 10 -4 is applied at the very 
first step of the run, which allows spiral modes to grow from the background noise as the 
disk cools. 
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3.2. The Evolution 

Figure 2 shows number density contours in the equatorial plane of the RC simulation. 
The RC disk undergoes the same four evolutionary phases described for CCT in Paper II, 
namely the axisymmetric, burst, adjustment, and asymptotic phases. In the axisymmetric 
phase, the disk shrinks slightly while cooling dominates over AV heating. The outer disk 
produces a dense, transient ring, and it becomes part of the spiral arms once instability 
sets in. The spiral structure develops between 2 and 3 ORPs of evolution, and the disk 
nearly doubles in size during the burst phase. The first expansion happens right at 4 ORPs, 
but, during adjustment, the disk clearly oscillates and reaches comparable sizes at 6, 7.7, 
and 10 ORPs. Subsequently, in the asymptotic phase, the disk maintains roughly the same 
outer radius. Often, one long spiral arm dominates the structure, and the disk has a more 
lopsided appearance than the CCT disk's analyzed in previous papers. Because the central 
star is kept fixed at the grid center, the dynamics of the these one-armed structures may 
not be accurately treated. We are developing routines to relax this constraint in the future 
by calculating the star's motion explicitly. Methods like moving the star to a location that 
brings the center of mass to the center of the grid (Boss 1998) are not employed because this 
might incorrectly treat the dynamics as well. Table 1 lists several properties of the initial 
and final disk at the equatorial plane for 10, 30, and 50 AU in radius. 

Table 1: Characteristics of the RC disk at and 16 ORPs. Here -Rdisk is the outer radius of 
the disk. 



t(ORP) 


-Rdisk 


Grid size (r, 0, z) 


r (AU) 


T(K) 


Q 


n(cm 3 ) 


E(g cm 2 ) 








10 


102.5 


7.4 


1.5(12) 


207.9 





40 


256,128,32 


30 


45.9 


1.6 


4.5(11) 


141.6 








50 
















10 


24.3 


1.8 


8.5(12) 


457.8 


16 


65 


512,128,32 


30 


8.4 


1.5 


5.2(11) 


49.8 








50 


3.4 


2.5 


4.8(10) 


11.9 



3.3. Mass and Density Distribution 

As expected, the structure of the disk changes significantly once the gravitational in- 
stabilities develop. Large amounts of mass move radially as the spiral pattern becomes 
nonlinear. We calculate average mass fluxes by differencing the total mass inside cylinders 
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Fig. 2. — Evolution of the RC disk. All the images show number densities at the 
equatorial plane on a logarithmic color scale that stretches four orders of magnitude 
with red representing the highest density and blue the lowest. The squares enclose the 
disk at maximum size, 170 AU on the side. A movie of the simulation is available at 
http://westworld.astro.indiana.edu/ under the Movies link. [Journal version: Evolution of 
the RC disk. All the images show number densities at the equatorial plane on a logarithmic 
gray scale that stretches four orders of magnitude with black representing the highest density 
and white the lowest. The squares enclose the disk at maximum size, 170 AU on the side. A 
movie of the simulation is available at http://westworld.astro.indiana.edu/ under the Movies 
link.] 

at two separate times and dividing by that time interval. This method should represent the 
average mass fluxes in the disk as calculated by the code's second-order flux scheme. Post- 
analysis methods for calculating the M's prove to be unreliable because the fluctuations in 
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the instantaneous mass fluxes are much larger by a factor of a few to about 10 than the net 
mass flux and are highly variable. As discussed in §5.1, the picture of mass slowly diffusing 
between different locations is incorrect for the RC disk. Between 2 and 4 ORPs, most of 
the net inward mass transport happens between about 9 and 27.7 AU, while net outward 
transport occurs outside 27.7 AU (Fig. 3). The accretion rates inside 7 AU are negligible. 
The peak transport rates are 1.75 x 10~ 5 and 3.25 x 10™ 5 M Q /yr inward and outwards, re- 
spectively. During the adjustment phase, the mass transport decreases by about an order of 
magnitude with typical rates of about 10 -6 M Q /yr and is mostly inward over 10 to 40 AU. 
After 10 ORPs, which marks the beginning of the asymptotic phase, the inflow rate peaks 
at 8.5 x 10 -7 M /yr, averages about 2 x 10 -7 M /yr, and is mostly inward over 10-25 AU. 

Figure 3 also illustrates how the r -0 5 surface density profile of the initial disk is lost, 
and at 4 ORPs the disk approximately achieves its final distribution with the profile obeying 
a Gaussian of the form £ = £ o 10~( r//r ' e - ) (dashed curve). The least-squares fit in log-linear 
space with the dummy variable x = r 2 yields r e = 46.7 AU when fit between r = [20, 60] 
AU. However, on the intervals r = [20,43] AU and r = [43,60] AU the surface density 
profile seems to follow two different power laws £ ~ r~ v . For the inner interval, v = —1.93, 
and for the outer interval, v = —5.97. The Spearman correlation coefficients are R = 
—0.992, —0.986, —0.985 for the exponential, the inner power law, and the outer power law, 
respectively. Michael et al. (2006, in preparation) explore the effects of initial conditions on 
these final surface density distributions. 

The disk inside 20 AU is not well-described by any simple monotonic function. There 
is a ring at 7.5 AU, and another seems to be forming at 10.5 AU, both resembling those 
that appeared in the CCT simulations. These rings are obvious in the cylindrical mass plot 
(Fig. 4). The broad peak seen around 15 AU in the cylindrical mass plot is not a true ring 
but represents mass concentrations in spirals that come together in the 10.5 AU ring. The 
mass of the 7.5 AU ring continuously increases and reaches 5 Jupiter masses (Mj) at the end 
of the simulation. The 10.5 AU ring has a mass of 9.0 Mj at 16 ORPs. See §5.1 for more 
details and a discussion about the causes and possible significance of these rings. 

3.4. Heating and Cooling 

The top panel of Figure 5 shows the total internal energy U contained in the disk as a 
function of time. The same plot also tallies total energy losses due to cooling in the optically 
thin regions (atmosphere) and in the interior of the disk as well as the total energy gained 
from heating by shocks. The overall behavior of the total internal energy of this disk is 
similar to the CCT disk's behavior. It must be noted that the — PV ■ v term is not included 
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Fig. 3. — Global mass redistribution. The top panel shows the average mass transport rate 
calculated by differences in the total mass fraction as a function of radius for three different 
times: between 2 and 4 ORPs, 4 and 10 ORPs, and 10 and 16 ORPs. The red [gray] curve is 
scaled to the left ordinate while the blue [starred] and green [dark] curves are scaled to the 
right ordinate. The bottom panel shows the surface density as a function of radius for the 
initial disk (dotted) and the final state (solid), which reflects the density profile during the 
asymptotic state. The dashed curve shows that the density profile for r > 20 AU follows a 
Gaussian distribution as described in the text. However, the surface density profile can also 
be broken down into two power laws (blue [dark] curves). 
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Fig. 4. — Rings in the inner disk. This plot shows the total mass per grid cylinder (1/6 AU 
in width) at various times during the evolution. The rings are labeled with their respective 
masses. 



in Figure 5. Therefore, just adding the energy gains and subtracting the energy losses due to 
A, T, and V • F from the initial value of the internal energy will not result in the final value 
of U. Moreover, due to coarse resolution, some of the r < 2/3 emission at the boundary of 
the atmosphere and disk interior is effectively included in the V ■ F curve. 

The evolution of the Toomre Q(r) of the disk is shown in the bottom panel of Figure 
5. As in Paper II, Q is calculated using the full vertical (top + bottom) surface density and 
both the sound speed and the angular speed evaluated at the midplane, where Q replaces 
the epicyclic frequency in equation (1) for a nearly Keplerian disk. All these quantities are 
azimuthally averaged before Q is computed. At 16 ORPs, the average value between about 
14 and 43 AU is 1.5, with a standard deviation a of 0.14. Together, the internal energy 
curve and the Toomre Q curves suggest that, after about 10 ORPs, RC is in an asymptotic- 
like phase as seen in the CCT disk with the disk remaining marginally unstable and shock 
heating roughly balancing cooling. The distinction asymptotic- like is made because RC's 
cooling times continuously adjust, and its evolution is different from the CCT disk's in that 
disk properties change noticeably even in this phase. We estimate the evolution time t ev for 
the disk in the asymptotic phase by t ev ~ U/U and by t cv ~ M in /M, where M in is the mass 
contained between 20 and 26 AU, ceteris paribus. Both methods yield t ev ~ few xlO 4 yr; 
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Fig. 5. — Total energies (top) and Toomre Q (bottom) as a function of time. The green 
curve in the top panel shows the total internal energy in ergs. The other curves show the 
total cumulative energy gains and losses due to the cooling and heating processes, also in 
ergs. Plotted separately are the contributions by shock heating, diffusion in the disk interior, 
and radiative cooling in the atmosphere. The bottom panel shows the Toomre Q every few 
ORPs. 



according to these estimates, the asymptotic phase is short-lived compared with typical total 
disk lifetimes (e.g., Hartmann 2005). However, effects like grain growth would significantly 
alter the calculation. Results presented by Cai et al. (2006) indicate that grain growth in a 
disk where the dust remains well mixed can lead to increased cooling times, which lead to 
weaker GIs and to a slower evolution. 

The average effective temperature (Fig. 6) over the 10-16 ORP time interval, computed 
by taking the fourth root of the azimuthally averaged T e g, seems to follow a temperature 
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profile T = T o 10~ r/re . For RC, r e = 44.8 AU when fit over the region r = (3,50] AU; this 
result is fairly insensitive to the limits chosen for the fit and has a Spearman correlation 
coefficient R = —0.988 in the log-linear space. However, if we choose to fit only the region 
r = (3, 20] AU, a power law T ~ r -o.59 seems ^ Q approximate the profile well with an 
R = —0.957. This value of the exponent is primarily due to the hot inner disk boundary 
between 2.3 and 5 AU. 




r(AU) 

Fig. 6. — Average effective temperature vs. radius time averaged over 10-16 ORPs. The 
profile that best fits most of the disk is an exponential, T = T o 10- r/re for r e = 44.8 AU, 
delineated by the red [gray] line. However, the r = (3, 20] AU region appears to obey a 
power law T ~ r ~ - 6 shown by the blue [dark] curve. 

The series of side-view maps shown in Figure 7 illustrate density and temperature at 
the end of the run. The left column shows a vertical cut through the disk, while the right 
column shows the azimuthal average. Generally speaking, the disk is dense in the inner parts 
and is more diffuse and slightly flared in the outer disk. Most of the disk interior, i.e., the 
region inside the disk photosphere, is contained inside about 25 AU, and its half thickness 
is only about 1 AU in the vertical direction, shrinking considerably from the interior of the 
initial disk with a radial extent of 37 AU and a vertical half-thickness of 3 AU. This is also 
the region where the temperatures remain higher than 10 K. The rest of the disk has cooled 
to temperatures of only a few Kelvin. However, there are streaks of higher temperatures in 
the outer disk created by shocks. Shock heating, in general, has the shortest times in the 
upper layers (Pickett et al. 2000) when compared to the local orbital time, and this heating 
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is probably due to shock bores, which are the nonlinear outcome of spiral shocks in disks 
(Boley & Durisen 2006). Typical cooling/heating times are a few to tens of ORPs, which 
indicates that the local cooling limiters defined in equation (14a) do not hinder the evolution 
of the disk in the asymptotic phase. 



Azimuthal Slice Azimuthal Average 




10 11 12 13 0.6 0.8 1.0 1.2 1.4 



Fig. 7. — Vertical structure of the disk at 16 ORPs for one slice of the grid (left) and the 
azimuthal average of all the slices (right). Number density is represented in the top panels 
and temperature in the bottom. 



4. RADIATIVE COOLING VS. CONSTANT COOLING TIME 

The RC disk undergoes gravitational instability much sooner than the t coo i = 2 ORP 
CCT disk in Paper II. During the axisymmetric phase, the typical cooling time of the op- 
tically thick (tr oss > 2/3) interior disk is shorter than 2 ORPs, so it cools to instability in 
about half the time. Short initial cooling times are due to strong departures from uniform 
radiative equilibrium in the initial model. As the simulation progresses and the disk adjusts, 
the typical cooling times increase, and, at the end, an entire column radiates its internal 
energy in a matter of a few to several tens of ORPs. We compute the azimuthally averaged 
column-wise cooling time by t coo i = (e co i) / (A), where e co i = J °° e dz is the vertically inte- 
grated internal energy of the half disk in a given column and A = J °° (V ■ F + A) dz. The 
brackets indicate azimuthally averaged quantities. As shown in Figure 8, the azimuthally 
averaged column-wise cooling times generally increase with time and with disk radius, which 
is consistent with the argument that t co0 \ ~ T~ 3 r for large r and t coo i ~ ^~ 3K pianck f° r sma U 
r by Durisen et al. (2006) (see also Rafikov 2005). Also shown in Figure 8 are curves for 
^cooi = 6/f2 and t coo i = Once the disk moves away from its initial state, the cooling 

times never drop below the Rice et al. (2005) fragmentation limit (t C ooi = 6/0) . We discuss 
the significance of the t coo i = 25/Q curve in §5.1.2. 

Figure 9 shows the final states of RC and CCT at 16 and 23.5 ORPs, respectively. We 
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Fig. 8. — Azimuthally averaged column-wise cooling times as a function of disk radius for a 
snap shot at 0.6 ORPs (squares) and for the temporal averages between 0.6 and 1.9 ORPs 
(crosses), 7 and 8.1 ORPs (triangles), and 12.2 and 13.2 ORPs (light curve). The heavy- 
curve profile represents the average cooling times between 10 and 16 ORPs as a function 
of r. For the profiles that continue past 40 AU, the cooling times continue to increase and 
become well over 100 ORPs. The profiles for 0.6 and 0.6-1.9 ORPs are radially smoothed. 
The dotted curves represent i cool = 6/fl (lower) and t coo i = 25/Q (upper). 

choose to compare the disks with each other at 16 and 23.5 ORPs because, although CCT 
has been evolved for a much longer time period, RC evolves much more quickly. The high 
temperatures in the upper layers of CCT are due to heat deposited by shocks, and, because 
the cooling time is constant everywhere, the upper layers cool more slowly there than they 
do in the RC simulation. This illustrates the importance of shock bores in disks. Boley & 
Durisen (2006) demonstrate that shock bores can deposit large amounts of energy into the 
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Fig. 9. — Comparison of the azimuthally averaged vertical structure of the RC and CCT 
disks at 16 and 23.5 ORPs, respectively. CCT's surface is much hotter than RC's because 
the cooling time scale is spatially uniform, while heating by shocks is highly nonuniform. 



disk atmosphere and its middle layers. These high-altitude temperatures should make the 
disk convectively stable (see §5.3); however, since the energy is deposited in layers that can 
cool efficiently, shock bores could effectively keep a disk cooler overall by depositing energy 
from the post-shock region of spiral waves into the high-altitude gaseous layers. 

The instabilities are stronger in CCT, with an average Toomre Q of 1.44 (er = 0.24) 
compared with 1.50 (a = 0.14) for RC, each measured between 14 and 43 AU. To measure 
the nonaxisymmetric structure in both disks, as done in Cai et al. (2006), we compute the 
time-averaged, integrated Fourier amplitudes (A m ), where we define 

Am= JpM. 

J po raraz 

Po is the axisymmetric component of the density and p m is the total Fourier amplitude of the 
cosm0 and the sinm0 density component (see Imamura et al. 2000). Summed over m = 2 
to 63, the integrated Fourier amplitude for t coo i = 2 ORP CCT is about 2.7 when averaged 
between 21.4 and 23.4 ORPs, whereas the amplitude is about 1.6 for RC when averaged 
between 14 and 16 ORPs, which is consistent with CCT having stronger GI activity than 
RC in the asymptotic phase. This is expected given the substantially longer column-wise 
cooling times in RC. 
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5. 



DISCUSSION 



5.1. Mass Transport and its Locality 



A principal question (Pringle 1981; Laughlin & Rozyczka 1996; Balbus & Papaloizou 
1999; Gammie 2001; Lodato & Rice 2004; Mejia et al. 2005) is the following: Can mass 
and angular momentum transport in Gl-active disks be modeled by a local a prescription, 
or do the long-range torques that GIs produce make such an approximation misleading? If 
GIs have a global effect on mass transport in the disk, substructure caused by spatially and 
temporally variable accretion rates, such as rings or density enhancements, can be missed 
by assuming that the disk evolves like an a-disk. In this section, we quantify the effects of 
GIs on mass transport in the RC disk and compare that behavior with CCT's to address 
the locality of GIs in disks. 



GIs redistribute mass and transport angular momentum over most of the RC disk chiefly 
through the torque exerted on the disk by the spiral waves. Peak mass transfer rates are 
larger than 10~ 5 M /yr during the initial burst and asymptote to a few xlO -7 M /yr 
at later times, comparable with those of the CCT disk. To investigate mass and angular 
momentum transport in RC further, we compute the torque in the disk. Following Lynden- 
Bell & Kalnajs (1972), one may write the torque C due to the outer disk on the disk inward 
of some distance x from the origin by integrating the gravitational stress tensor T grav and 
the Reynolds stress tensor y Reyn over the surface of a cylinder; 



For the gravitational contribution, instead of using the stress tensor directly, one may inte- 
grate over the total volume by 



where $ is the gravitational potential. Furthermore, we are only interested in the z direction 
of the torque C z given by 



5.1.1. Torques and Modes 




(16) 




(17) 




(18) 



The contribution to the z-torque from the Reynolds stress is calculated by 




(19) 
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where Sv^ is the fluctuating component of the azimuthal velocity and 5v r is the fluctuating 
component of the radial velocity. The fluctuating components are calculated by subtracting 
the mean velocity, found by azimuthally and vertically averaging the velocities with a density 
weighting at each r, from the velocity in each cell for the corresponding r. A more detailed 
discussion of the analysis associated with equations (17) to (19) will be given in Michael et 
al. (in preparation). 

Figure 10 shows the time-averaged gravitational torque (heavy black) and gravitational 
plus Reynolds torque (light black) over 10-16 ORPs; the red curve is explained below. Su- 
perposed onto the plot is the average M(r) over the same interval. The ring at 10 AU 
(see Fig. 4) is formed where the total torque precipitously declines. This also agrees with 
the sudden drop in accretion rate, as discussed in §3.3, near the same radius, and indicates 
that the formation of this ring can be understood as the place where mass piles up due to 
a diminution of GI torques at this radius. It is unclear whether this ring would survive if 
our code included additional accretion mechanisms in the inner disk, e.g., an a viscosity 
caused by the magnetorotational instability (Balbus & Hawley 1991; Gammie 1996). This 
is a subject for future study. 

To complement this figure, we present a periodogram (Fig. 11) for the m = 2 Fourier 
component of the nonaxisymmetric density structure in RC as done in §3.3.2 of Paper II 
for CCT. The periodogram detects Fourier component m-armed density fluctuations that 
have coherent pattern frequencies over a large range of radii. The contours show power, 
where dark red is the strongest and light blue is the weakest, and the dark curves indicate 
corotation (Q = f2 pa ttcrn) with the gas and the inner and outer Lindblad resonances (fi = 
Q±K e /m). The resonance curves are noisy in the ring radii range, even after some smoothing, 
because the rings create strong pressure gradients, which affect the epicyclic frequency. There 
is significant power at corotation at a period of about 4/7 ORP near r = 26 AU. This 
corresponds to the mass transport inflow/outflow boundary, which also corresponds to the 
global maximum in the gravitational torque and a local maximum in the total torque in 
Figure 10. Moreover, there are swathes of power near periods of 1/2 and 1 ORP. The 
1/2 ORP swath has its outer Lindblad resonance at the inflow/outflow boundary, and the 1 
ORP swath has its inner Lindblad resonance at the inflow/outflow boundary. The alignment 
of the three major swathes of power with corotation, inner Lindblad, and outer Lindblad 
resonances, combined with coincidences between corotations and significant M's, as shown 
in Figure 11, suggests dominance of low-order modes in the mass and angular momentum 
transport. In addition, there is power on the corotation curve at 1/9, 1/6, and 1/4 ORP; 
the rings are between these locations. As suggested by Durisen et al. (2005), resonances may 
play a role in forming and maintaining the rings. 




Fig. 10. — Time-averaged gravitational and Reynolds torque profile (thin solid) and just 
the gravitational torque (heavy solid) for the 10-16 ORPs evolution time interval. We use 
the sign convention that a positive value indicates outward angular momentum transport. 
Superposed on the plot, unsealed to the ordinate, is the time-averaged mass accretion rate 
(dashed line). The major break between inflow and outflow occurs near the maximum of 
the gravitational torque (r = 26 AU), but the highest mass inflow occurs inside the global 
maximum of the total torque at r = 16 AU. The red [gray] curve is the torque profile derived 
from the M's for the 10-16 ORP interval. See section 55.1.3 for more detail. 



The three major power swaths are also noticeable in the periodograms for m = 3, 4, and 
5, but they become less prominent with increasing m. Moreover, the periodograms become 
more noisy with increasing m in the sense that there are many other swaths of power present. 
Therefore, we present in Figure 12 an amplitude spectrum for the A m 's defined in equation 
(15) for m = [1,63]. Clearly, the low-order m's have the largest amplitudes. The m = 1 
amplitude is probably larger than what it would be if we would allow the star respond to 
the background potential, and the extra power in the m = 1 mode probably affects the 
disk dynamics in some respects but not in the gross behavior, which we believe based on 
the consistent results between the power in the m = 2 mode, the periodograms, and the 
mass transport. For m = [2, 63], the data are fit by the model curve A m ~ (m 2 + m 2 ,) -1,64 , 
where m = 7.46. This curve highlights that A m asymptotes to a power law for large m 
values. Along with the noise in the periodograms, the A m profile suggests that GIs lead 
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to a cascade of power among many modes resulting in gravitoturbulence (Gammie 2001). 
In fact, the radial wavenumber power spectrum shown by Gammie (2001), which excludes 
low-order modes because a shearing box approximation is used, appears to be consistent 
with our power spectrum for m > 10. Understanding the slope of this spectrum for large m 
values and determining the uniqueness of this profile are topics for future study. 



5.1.2. Is Angular Momentum and Energy Dissipation Accurately Described by an a 

Prescription ? 

The torque can also be used to calculate a Shakura & Sunyaev (1973) a for the disk. 
Define 

a = dlnQ ~' (G r<p ) + (R r4> ) ^ 
dlnr (^ c s) 

where (G r( p) is the average vertically integrated r-0 component of the gravitational stress 
tensor, (R r <f,) is the vertically integrated Reynolds stress tensor, £ is the surface density, c s 
is the midplane sound speed, and the brackets indicate that the quantities are azimuthally 
averaged (see Gammie 2001; Lodato & Rice 2004). We may use the torque C z to find the 
stress from equations (16), (18) and (19), 

C z = 2nr 2 ((G r ^ + (R^}). (21) 

Using equations (18)-(21), we plot the time-averaged effective a for the 10-16 ORP time 
period in Figure 13. On the same figure, we show the a that is expected for a disk when 
entropy generation by gravito- and hydrodynamic turbulence provides enough heating to 
exactly balance cooling, namely, the a given by equation (2). For the stretch between about 
20 and 35 AU, Figure 8 shows that the t coo i = 25/S1 curve roughly overlaps the average t coo i 
for the 10-16 ORP time interval. If we take t C ooi^ ~ 25 over that interval and calculate the 
a expected from equation (2), we find log a = —1.9 for the strongly self-gravitating case 
and logo? = —1.6 for the negligible self-gravity case (see text below and equation [2]). We 
also calculate an a based on the bold t coo i curve in Figure 8 with equation (2), which is 
averaged over 10-16 ORPs, and by assuming the negligible self-gravity case. In Figure 13, 
the a expected from equation (2) with the t coo i curve (heavy dashed curve) is consistent with 
the a derived from the gravitational torque with equation (20) (heavy solid curve) between 
15 and 60 AU in r. This a is of the same magnitude as the a used in SED fitting for real 
disks (e.g., Hartmann et al. 1998), but only in the 20-35 AU region is the effective a, which 
is of the order 10~ 2 , roughly constant. 
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Fig. 11. — Periodogram for the m = 2 mode. Even through the noise, three large swathes 
of power are noticeable at about 1/2, 4/7, and 1 ORP, which cross the outer Lindblad, 
corotation, and inner Lindblad resonances, respectively, at the mass inflow/outflow bound- 
ary. Moreover, there is significant power at about 1/9, 1/6, and 1/4 ORP at corotation. 
The radii corresponding to these local power maxima straddle the rings. The dotted line 
indicates where the mass accretion rate almost returns to zero and corresponds to the coro- 
tation resonance for the swath of power near 1/2 ORP. This line also represents where the 
exponential surface density profile abruptly ends with a sudden plateau in the overall profile 
and the formation of rings (Fig. 3). 

5.1.3. Are the Torque Profiles Consistent with the Reported M's? 

To check the consistency of the reported M's with the torques in Figure 10, we calculate 
the torque that is required to produce the M's shown in Figure 3 for the 10-16 ORP interval. 
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r. 



With azimuthally averaged quantities so that each quantity only depends on 

^ = ^< c »>' (22 > 

Here, we have assumed that dlnf2/dlnr = —3/2, which is a very good approximation for 
most of our disk. Equation (22) simply reflects the drift of mass due to changes in its orbital 
radius caused by torques. By integrating this equation, we calculate the total torque required 
to produce the M's measured in the RC disk. 

We can view the consistency check from a different vantage point by using the torque 
from equation (22) to derive a corresponding a. We relate the torque and a by letting 
v = ac s h, by assuming that the disk scale height h ~ c s /Q, where the midplane values are 
used for c s and Q, and by using equations (20) and (21) in equation (22), which yields 

(23) 



The resulting torque profile from equation (22) and the a profile from equation (23) 
are shown as red curves in Figures 10 and 13, respectively. This a profile is reasonably 
consistent with the a profile derived from the gravitational torque alone for the 5 to 50 AU 
internal. Over that interval, the M-derived a's are slightly larger than the gravitational 
torque-derived a's; this is probably due to ignoring constants of order unity in some of the 
relations. For r > 50 AU, the a's deviate strongly. 

In Paper II, we reported that the effective a's seen in our disks were an order of magni- 
tude higher than those given by equation (2). This is probably due to use of the mass flux in 
a steady-state, namely M = 37rz/£, to compute the effective a's from the M's. We suspect 
that the anomalously high effective a's in Paper II are incorrect for this reason; Michael et 
al. (in preparation) will examine the effective a's in constant t coo i disks with equations (20) 
and (23). 



5.1.4- Summary 

The M's in the RC disk are primarily produced by global modes, as shown by Figures 
10, 11, and 12, but the a profile derived from the torques and mass transport is consistent 
with a locally applicable a(r,t). This is likely due to the close placement of corotation 
radii for three strong m = 2 modes over the 20-35 AU region. As predicted by Balbus & 
Papaloizou (1999), a disk that is dominated by global modes will behave like an a-disk as 
regards angular momentum transport and energy dissipation near the corotation radii of 
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low-order modes. We note that we are unable at this time to determine whether the energy 
dissipation in RC is consistent with some locally defined a(r,t), but we do find that the 
angular momentum transport is consistent. We are currently addressing this issue, and plan 
to discuss it in a future paper. Moreover, the disk deviates from an a-disk behavior for 
r > 50 AU, which is far from the main m = 2 corotation radii; this is consistent with Balbus 
& Papaloizou (1999). 

Based on the the torque profile, the a profile, the t coo i^ profile, and the Fourier amplitude 
spectrum of the RC disk, we conclude that a local description for mass transport is consistent 
with an a-disk near the corotation radii of low-order modes. We stress, however, that the 
evolution is not described by a constant a. Using an a(r, t) prescription is only possible if one 
knows the correct t C ooi' s f° r the asymptotic phase, which requires hydrodynamics calculations 
(see Johnson & Gammie 2003), and if one can estimate the location of corotation radii for 
low-order modes a priori. In any case, the low-order modes lead to significant radial mass 
concentrations, e.g., Figure 4, which would be lost in a smooth a prescription. 

5.2. Fragmentation and Rings 

The RC disk develops a series of rings as massive as the ones seen in the CCT disk. 
Matter is gathered at these radii and the rings are fed continuously, as seen in Figure 4. 
This, along with the discussion in §5.1, implies that ring formation may be a robust result 
as long as there is a region of the disk with negligible torque or a region with kinks in the 
torque profile. 

Recently, while conducting tests of this radiative scheme against a new scheme that 
implements rays for part of the radiative solution, it was noticed that the portion of the disk 
inside the first ring may be staying hot for numerical reasons. This does not affect points 
made here about the rings, i.e., rings can build up where the torques precipitously decline. 
In this case, the torques probably decline because, inside 7 AU, the disk is stable to GIs 
due to numerical heating. However, there are several physical heating mechanisms that are 
excluded from this calculation, and we expect to see the same results if the inner 7 AU of 
this disk were to be kept hot for a physical reason. We are currently running a series of 
calculations that do not have such a large radial range so that we may avoid poor numerics 
in the inner portion of the disk. 

Durisen et al. (2005) speculate that concentrations of mass may contribute to planet 
formation by accelerating core accretion (see also Pickett & Lim 2004, Paper II §4.3). This 
might be the only way that GIs in RC could produce planets. Typical cooling times in RC 
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Fig. 12. — The A m 's vs. m in log-log space. The squares represent the data derived from 
the RC disk and the heavy line is the best fit curve through those data, excluding m — 1. 
The error bars represent the rms of the residuals and indicate the degree of fluctuation 
in the spectrum. There is a noticeable wiggle in the linear portion of the curve, so the GI 
amplitude spectrum is probably more complex than described here. This is a topic for future 
investigation. 

increase from fractions of an ORP to a few and even tens of ORPs with time (see Fig. 8). 
Initially, t C ooi^ J$ 6 for several radii, but quickly (within 1 ORP) the cooling times increase, 
and the disk becomes strongly stable against fragmentation at all radii, because, as shown 
in Figure 8, t C ooi^ ^> 6 (Gammie 2001; Johnson & Gammie 2003; Rice et al. 2005). This 
is also consistent with the stress fragmentation criterion (Rice et al. 2005), which requires 
that the effective a in nonfragmenting disks be < 0.06. The Toomre Q drops below 1.0 
only just before the burst phase at about r = 30 AU. To ensure that no fragmentation was 
overlooked due to resolution effects, the burst phase was rerun with 4 times as many cells in 
the azimuthal direction (512 instead of 128). No significant structural change was observed. 
Dense clumps do form during the burst phase at the intersections of spiral structures, but 
they last only a fraction of an ORP as in Papers I and II. In addition to the burst phase, the 
asymptotic phase was extended for one ORP with the same high resolution. Qualitatively, 
the disk structure remains the same and no signs of fragmentation are detected. 
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Fig. 13. — Effective a based on the total torque profile in Figure 10 (thin solid curve) and the 
gravitational torque only (heavy solid curve). The red [gray] curve indicates the a profile one 
would need to yield the M's for the 10-16 ORP interval shown in Figure 3. The dash-dot lines 
are effective a's predicted by Gammie (2001) for t C ooi^ = 25 for the strongly self-gravitating 
case (lower) and the nonself-gravitating case (upper). The heavy dashed curve indicates the 
predicted a profile based on the t coo i's from the 10-16 ORP interval, which are shown in 
Figure 8, with the assumption of negligible self-gravity. 

The long cooling times in RC are consistent with the predictions of Rafikov (2005), with 
cooling times in the Nelson et al. (2000) simulations (see §4.2 Durisen et al. 2006), and with 
the cooling times due to radiative cooling alone in the Boss (2001, 2002, 2005) simulations. 
However, Boss also claims to see convection in his simulations, which decreases the cooling 
times in his disks enough to allow for disk fragmentation and the formation of multi- Jupiter 
mass clumps. We see no signs of convection in RC during the Gl-active phases. 

5.3. Convection 

Lin & Papaloizou (1980) (see also Ruden & Pollack 1991) showed that a vertically con- 
tracting but otherwise quiescent protoplanetary disk will be convectively unstable when the 
optical depths are large, in the Rosseland mean sense, and when (3 > (37 — 4) / (7 — 1), 
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where (3 is defined as the exponent of the temperature dependence of mass absorption co- 
efficient Xross, i-e., Xross = XoT 13 - For a 7 = 5/3 gas, one should expect convection when 
(3 > 1.5, and for 7 = 7/5, one should expect convection when (3 > 0.5. In our calculation, we 
find that « 2 for T « 30 — 100 K, which is the temperature range in most of the interior 
(see Appendix A, Fig. 16). Indeed, we see convection in the axisymmetric phase of the disk 
evolution during which the disk is, for the most part, undergoing quasistatic gravitational 
contraction in the vertical direction. 

We search for convection in the simulation by looking for convective "cells" that are 
associated with regions where the entropy gradient is negative, i.e., <9(P/p 7 ) jdz < 0. We 
use the entropy gradient instead of the Schwarzschild criterion for convective instability, 
namely V = <91nT/<91nP > 0.4 = V a( j for a 7 = 5/3 gas, because the Schwarzschild 
criterion assumes a vertically hydrostatic background, which is not guaranteed in a dynamic 
disk. 

Figure 14 shows superadiabatic regions in the disk at 1.2 ORPs (top panel) and at 
10.6 ORPs (bottom panel). Both panels show heating due to artificial viscosity by gray- 
filled contours, velocity vectors scaled to each axis but with the heads on the left panel being 
slightly larger for clarity, and heavy curves representing density contours. The superadiabatic 
regions in the right panel are delineated by blue thin curves. The left panel omits these 
contours to show the velocity vectors more clearly; almost the entire region inside the inner 
most density contour is superadiabatic. In the left panel, the average Mach speed in the 
convective eddies (v z /c s ) ~ 0.06, where (v z /c s ) = f p \v z \ jc s dV/ f p dV and the integrals 
are evaluated over the volume spanning between 15 and 25 AU in r. The value of (v z /c s ) 
fluctuates between 0.01 and 0.15 when each annulus of grid cells is evaluated separately (Ar = 
1/6 AU), and some of the convective eddies result in nontrivial compressional heating through 
artificial viscosity. To evaluate energy transport by convection, we estimate the convective 
flux F c = — l/2c p pTv z £ (dlnT/dz — VaddlnP/dz) (Lin & Papaloizou 1980) through each 
cell in the volume described above. Here, the mixing length I = min (z, P/VL\zp) and c p 
is the specific heat at constant pressure. By dividing the total internal energy within the 
half disk between 15 and 25 AU in r by the total convective energy loss rate for that same 
region, we find that the convective cooling time is about 1 ORP and that it is comparable 
to the radiative cooling time. However, this method measures the convective flux based 
on a superadiabaticity estimate according to the formalism of Lin & Papaloizou (1980). 
For a second measurement, we calculate the energy carried by convective motions through 
a plane that cuts through the volume at z « 1 AU. The convective flux through the ith 
cell is F c |j = pv z c p AT\i, where AT is the difference between the actual temperature at 
the cell center and the azimuthally averaged temperature for that r and z. We find a 
convective cooling time of about 2 ORPs with this method. According to either estimate 
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of the convective flux, both of which are crude and uncertain, convection is efficient at 
redistributing energy in the disk's interior. However, these high convective fluxes may be due 
to a combination of the random perturbation to the initial constant vertical entropy profile, 
of our seeding superadiabatic regions at the interior/atmosphere interface (see Appendix 
B), and of our exclusion of irradiation, which tends to produce stabilizing temperature 
stratification. Moeover, this energy ultimately must be radiated away, so the cooling times 
in Figure 8, based on energy actually radiated, accurately reflect disk cooling times. 

The convective cells that are established in the axisymmetric phase are disrupted during 
the burst phase, and, for the subsequent Gl-active phases, we do not see convection. We 
believe this is due to two effects: (1) The dynamics in the disk, namely spiral waves, disrupt 
any establishment of convective cells due to the large fluctuations in mass transport. (2) 
These spiral shocks lead to shock bores, which in turn lead to high temperatures at moderate 
and high disk altitudes (Boley & Durisen 2006, see also Fig. 7, here). These high-altitude 
shocks should make the disk convectively stable by keeping the entropy gradient of the 
disk positive or essentially zero in the vertical direction. The right panel in Figure 14 
shows this well. Despite the presence of superadiabatic gradients (see appendix B for an 
explanation why there are strong superadiabatic regions at mid-disk altitudes), convective 
cells are absent. The only vertical motions associated with the flow are due to shock bores 
along the spiral shocks or other waves. 

Because shock bores deposit energy into the upper layers of the disk, it is reasonable to 
ask whether shock bores themselves provide convective-like cooling. Although this point is 
still unclear, we believe not. Shock bores effectively reduce the temperature in the post-shock 
region as the jumping gas expands vertically. This removes thermal energy in the post-shock 
region and deposits that energy in the upper layers via waves. Therefore, we do expect shock 
bores and possibly other wave effects (e.g., see Lubow & Pringle 1993; Lubow & Ogilvie 1998; 
Ogilvie 1998) to enhance disk cooling when large-scale spiral shocks are present in the disk. 
Although energy may be effectively removed from the post-shock region, the interior of the 
disk can still only cool by radiative diffusion because of the positive vertical entropy gradient 
established by the shock bores. From this point of view, shock bores and other wave effects 
limit the efficacy of spiral waves in heating the disk. Distinguishing between shock bores and 
convection is not a point of semantics. Shock bores are born of large-scale shocks in a disk, 
while thermal convection and the criterion for convective instability are usually described in 
the context of a disk in vertical hydrostatic equilibrium. 
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Fig. 14. — The dark lines indicate density contours, the gray shading indicates heating from 
bulk viscosity, and the thin blue/purple lines indicate superadiabatic regions, which are 
defined wherever the vertical entropy gradient is negative. The arrows show the direction of 
gas flow in this azimuthal slice. Top: An r-z slice through the disk during the axisymmetric 
phase (1.2 ORPs) when convection occurs. The entire region below the innermost density 
contour is superadiabatic, so the contours are omitted. Bottom: A slice through the disk but 
at 10.6 ORPs and for a larger region. Convection is absent in the superadiabatic regions; gas 
that is rapidly moving upward is related to shocks in the disk and are partly shock bores. 
The velocity vectors have a different scaling because the velocities are large. 



5.4. Spectral Energy Distributions 

As emphasized by Nelson et al. (2000), an important test of physical relevance is a 
comparison between SEDs of observed systems and the derived SED from a simulation. 
We attempt to construct what the RC disk would look like to an observer at some large 
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distance d looking at the disk face-on. For simplicity, we assume that each (r, <fi) column 
of the disk emits radiation according to a black body law at its effective temperature over 
its surface area. Because we use a cylindrical grid, the ith area element has a solid angle 
dQi = Tidrdcf)/ d 2 . The specific flux, or flux density, then can be tallied by 



where B v is the Planck function and Tj is the effective temperature of the ith cell. To avoid 
using distance, we choose to express the SED in terms of ^d 2 v¥ v . As a basis for comparison, 
we adopt the approach of Nelson et al. (2000) and define our fiducial to be an SED derived 
from a disk with a temperature law r -0,6 , which is the median best fit law to the T Tauri 
disk sample presented in Beckwith et al. (1990). Figure 15 shows the SEDs derived from 
RC for three different stages in its evolution. The short dashed line delineates the SED for 
the disk near its brightest period during the burst phase (2.4 ORPs), when the luminosity 
is nearly 19 L & as integrated between 10 10 and 10 15 GHz. The long dashed line delineates 
the SED for the disk as it enters the asymptotic phase (10 ORPs), and the solid black line 
indicates the SED for the disk at the end of the calculation (16 ORPs). The star, which 
is assumed to have an R = 2R Q and a T e fj = 4, 000 K, is included in the SED profile; the 
SED has dips in specific luminosity as it nears the star because it has a 2.3 AU hole and is 
missing a contribution from an inner and hotter portion of the disk. The red lines indicate 
fiducial SEDs based on assuming a T c g ~ r -°- 6 temperature profile for a 0.0033, 0.01, and 
0.033 L & disk and integrating between and 60 AU. The fiducial SEDs also have a slight dip 
in their profile just before transitioning to the stellar portion of the SED due to discretizing 
the temperature profile into grid cells 1/6 AU wide in r. The blue line is a fiducial SED 
calculated for a luminosity of 0.0024 L & , which is the disk luminosity at 16 ORPs, in the 
same way as the other fiducials but with a 2.3 AU hole. Even though the actual effective 
temperature profile for RC is an exponential (see Fig. 6), the SED for RC is observationally 
consistent with a profile T e fj ~ r -o-6 Decause the inner 20 AU of the disk closely follows a 
T c s ~ r -°- 59 . Although GIs appear to lead to an exponential T e g profile for this disk, stellar 
irradiation (D'Alessio et al. 2001) probably keeps the outer regions of the disk warmer than 
what is modeled here inasmuch as that region of the disk is flared (see Fig. 9). 

The SEDs of RC also suggest that an FU Ori outburst may be the signature of a disk 
becoming gravitationally unstable. The RC disk reaches a peak luminosity of 19 L & at about 
2.4 ORPs, which makes RC approximately 10,000 times brighter in the burst phase than in 
the asymptotic phase. Furthermore, the disk has a luminosity of only about 0.01 L Q around 
2.1 ORPs, so the disk luminosity increases by a factor of about 1,000 in 80 years. The 
increase in luminosity is sudden and may be much shorter than we are reporting because 
our time resolution for determining disk luminosity is only about 0.3 ORPs due to data 
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Fig. 15. — The solid line indicates the SED for RC at the end of the calculation (16 ORPs). 
The short dashed and long dashed curves delineate the SEDs for RC at 2.5 ORPs and at 10 
ORPs. The red [gray solid] curves represent fiducial SEDs for 0.0033, 0.01, 0.033 L Q for a 
disk that extends from 1/6 AU to 60 AU and has a T e g ~ r~ q for q = 0.6. The blue [dark] 
dashed curve indicates the SED for a disk that is truncated at 2.3 AU with a q = 0.6. 

handling restrictions, which also means we may be missing the peak brightness. GIs occur 
on a dynamical time scale and the steeping of spiral waves into shocks may be rapid enough 
to be consistent with observed FU Ori outbursts, which have luminosity rise times between 
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a few to tens of years (Hartmann & Kenyon 1996). In addition, the decay time scale for RC 
is about one ORP (about 250 yr), which is consistent with observed decay times. 

Although the burst phase is a promising explanation for FU Ori outbursts, one may 
certainly ask the question: Is the burst phase an artifact of our initial conditions? The answer 
is unclear. The burst phase is the transition from a globally stable disk to a globally unstable 
disk. We speculate that the burst we see in this disk might be applicable to a real disk and to 
FU Ori outbursts under the following scenario: A real disk might form a dead zone (Gammie 
1996; Armitage et al. 2001) wherein GIs may undergo episodic bursts. These bursts of GIs in 
an annulus can create strong shocks over a large Ar (Boley & Durisen 2005) and produce FU 
Ori outbursts with very short rise times. The location of the unstable annulus determines 
the rise time for the outburst. Once the disk becomes unstable, the shocks extend deep 
into the inner disk, which is excluded from this calculation and should lead to even brighter 
bursts with more luminosity at shorter wavelengths. The global torques exerted on the disk 
as a result of the burst are about 10 times greater than in the asymptotic phase, which 
results in effective a's that exceed lO^ 1 without leading to long-term fragmentation. The 
strong shocks and redistribution of mass lead to the restabilization of the unstable annulus 
and of the disk as a whole, and the increase in the envelope irradiation will also work to 
stabilize the disk (Cai et al. 2006). Furthermore, an annulus of material with a mass around 
10 Jupiter masses is expected to be completely redistributed in the nebula in a few hundred 
to a few thousand years, assuming a mass flux of between 10 -5 and 10 -4 M yr -1 . More 
detailed modeling is necessary, but we believe that there is enough evidence to suggest that 
FU Ori outbursts may be the signature of the sudden onset of GIs in protoplanetary disks. If 
this is true, then FU Ori outbursts may also signal the annealing of dust and the formation 
of chondrules in protoplanetary disks by spiral shocks (Wood 1996; Boss & Durisen 2005; 
Boley et al. 2005; Boley & Durisen 2005). 



5.5. Fate of the Disk 

The internal energy profile in Figure 5 shows that, after the burst, the disk transitions 
to a self-sustaining marginally unstable disk, but, as discussed in §3.4, that phase can only 
be maintained for a few xlO 4 yr. Because this is an incomplete model, i.e., we exclude the 
inner disk, irradiation, possible infall, dust settling, grain growth, and other mass transport 
mechanisms, it is difficult to comment on the disk's ultimate fate. Nevertheless, we speculate 
about the bigger picture, and hope that it might serve to link several areas of disk research. 

The time scale for grain growth from a maximum size of about 1 /im to about 1 mm 
is expected to be between the orders 10 3 -10 5 yr (Haghighipour 2005). This time scale is 
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commensurate with the lifetime of the asymptotic phase of the RC disk, ceteris paribus, 
and it therefore must have an effect on the disk's evolution. Even though our disk has a 
relatively high effective a, it is unclear whether there is enough turbulence in the disk to 
keep the larger mm-size grains from settling to the midplane. To proceed, we assume that 
no significant settling takes place while the GIs are as active as reported in §5. Under this 
scenario, the larger grain size will increase the opacity in cooler portions of the disk where the 
typical photon wavelength is long, while the inner disk will see a decrease in opacity where the 
typical photon wavelength becomes shorter than the maximum grain size (see, e.g., D'Alessio 
et al. 2001). We expect that, in the cool parts of the disk, the cooling times will increase 
as the maximum grain size grows; this is consistent with the simulations presented by Cai 
et al. (2006). Longer cooling times lead to weaker GIs and to slower evolution. As grains 
grow, the GIs may become so weak that they no longer provide enough turbulence to counter 
settling. Because the dust settling time is similar to the grain growth time (Weidenschilling 
1997; Haghighipour 2005), dust settling is rapid. As the dust settles to the midplane, the 
opacity for most of the disk decreases, which would decrease the cooling times. This would 
lead the gas disk toward stronger instability. If the instability is strong enough, the disk 
may undergo multiple burst-like phases, during which the effective a's ~ 10 _1 . Because 
rings are a pressure maximum that lead to enhanced dust to gas ratios (Weidenschilling 
1977, 1997; Haghighipour & Boss 2003; Rice et al. 2004; Durisen et al. 2005; Haghighipour 
2005), if rings form at the outer radius of a dead zone, one might expect grain growth 
and dust settling to happen there first. Whether the entire outer disk becomes unstable or 
just a large annulus, these bursts could re-mix the grains vertically and radially in the disk 
(Boley et al. 2005; Boley & Durisen 2006) and lead to collisional grain destruction and to 
the processing/reprocessing of crystalline material. Because bursts redistribute the mass in 
the disk, the lifetime of the disk increases. 

An opacity increase from grain growth can lead to a cessation of GI activity. The 
cessation of strong gravitoturbulence permits dust settling if there are no other means for 
producing strong turbulence. This dust settling might decrease the opacity enough that the 
cooling times become short enough to induce another disk burst. This sequence may be 
able to explain FU Ori outbursts and their duty cycle, drive shocks for chondrule formation 
(e.g., Desch & Connolly 2002, for a discussion on chondrule formation via shocks), and might 
extend the RC disk lifetime to ~ 10 6 yr. Although speculative, we believe this provides an 
avenue for further study. 
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6. CONCLUSIONS 

Our simulations of a protoplanetary disk show that (1) radiative cooling times are 
too long to permit disk fragmentation for a 7 = 5/3 when realistic opacities are used (2) 
mass and angular momentum transport are dominated by low-order modes, (3) angular 
momentum transport in gravitationally unstable disks with realistic cooling is consistent 
with the transport expected for an a-disk over a nontrivial region of the disk, (4) evolving a 
Gl-active disk as if it were an a-disk may be impossible to implement correctly because it 
requires knowing the tool's and the distribution of corotation radii for the low-order modes 
a priori, (5) the GI activity in our disk is weaker than in the constant cooling time disks 
of Mejia et al. (2005) because the cooling times are typically longer, (6) the duration of 
the Gl-active phase in this disk is expected to continue for a few xlO 4 yr, (7) convection 
is present during the axisymmetric phase but convective cells are completely disrupted by 
Gl-activity, and (8) features in the torque distribution can lead to the development of rings 
near the edge of a Gl-active region. We also find that our simulation yields an SED that is 
compatible with observed SEDs and that the GI burst phase may correspond to an FU Ori 
outburst-like event. 

The behavior of gravitational instabilities in disks is extremely sensitive to the handling 
of thermal physics, and researchers must use extreme caution when treating radiative bound- 
ary conditions. We find most of the disk volume is in the atmosphere, which is optically thin 
to its own radiation, as also found in the Cai et al. (2006) simulations. Therefore, proper 
disk evolution depends on the treatment of radiative cooling at low optical depths. Cooling 
times tend to increase past the initial phase to the point that a vertical column can take 
from several hundreds to tens of thousands of years to radiate its internal energy. These 
time scales agree with those seen by Boss (2002) due to radiative cooling alone and are too 
long to produce planet formation by disk fragmentation. Convection could provide shorter 
cooling times (Boss 2002), but convective motions are absent in the simulations presented 
in this paper and in Cai et al. (2006) during Gl-active phases. Even if the shock bores that 
produce vertical motion in our disk enhance the cooling, the cooling times are still long. 

Assumptions about the physical properties of dust grains that are mixed with the gas 
in these disks play an important role in determining the cooling rates of the entire disk. 
The problems of dust settling (Calvet et al. 2005) or aggregation of solids (Rice et al. 2004) 
are not addressed in these models, but they will be significant factors to consider in future 
simulations. Moreover, a 7 = 7/5 gas disk may behave differently from what is modeled 
here (Rice et al. 2005). Including all the details necessary to handle the flow of energy in 
circumstellar disks not only requires a large amount of future human effort and computing 
time, but it is also requires better constraints on the physical parameters of real disks. We 
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are now developing a new radiative scheme that couples ray tracing with radiative diffusion 
in an effort to improve coupling of optically thick and thin regions. 
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A. OPACITIES 

All the opacities used in the radiative cooling simulations were obtained from Paola 
D'Alessio. The opacities are similar to those by Pollack et al. (1994), but with a few modifi- 
cations (see D'Alessio et al. 2001 for details). The main contributors are grains of H2O ice, 
silicate (olivine [Fe,Mg] 2 Si0 4 and orthopyroxene [Fe,Mg]Si0 3 ), organic (containing C, H, O, 
and N), and troilite (FeS) (Figure 16). The dust grains are assumed to be spherical, and 
their size distribution is a power law 

dn = noa~ s da, (Al) 

where a is grain radius, no a normalization constant, and s a free parameter. The 
opacities used in all the simulations assume s = 3.5, which best fits observed ISM extinctions 
for grains of various compositions (Mathis et al. 1977; Pollack et al. 1994). The minimum and 
maximum grain radii are 0.005 and 1 urn, respectively, although opacity tables for different 
maximum grain radii, a max , are available. While these grain sizes are more representative 
of interstellar dust (Draine & Lee 1984), they have been widely used to model composition, 
abundances, and physical properties of the Solar Nebula and circumstellar disks in general 
(e.g., Pollack et al. 1985, 1994; Calvet et al. 1991; Henning & Stognienko 1996; Bell 1999; 
D'Alessio et al. 1999). This appendix gives the definitions of the various mean opacities 
used in the 3-D hydrodynamics code, summarizes their properties, and briefly compares how 
different grain sizes can affect disk structure and energy transport. For details on the specific 
opacities that are used to calculated the mean opacities described in this appendix, we refer 
the reader to D'Alessio et al. (2001). 
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Since the 3-D hydrodynamics code does not compute frequencies or wavelengths, the 
opacities used in the calculations of the heating and cooling terms in (§2.1) are the Planck 
and the Rosseland mean opacities integrated over all frequencies, defined as 



and 



^Planck 



J °° K V B U &V 
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respectively, where k v represents only true absorption while \v constitutes total extinction, 
i.e., absorption plus scattering and where B v is the Planck function. Rosseland opacities 
are used in optically thick regions, where the diffusion approximation is valid and the mean 
free path of a photon is much smaller than the thickness of the gas. Planck opacities better 
characterize regions where photons are likely to escape without interacting with the medium 
(e.g., Mihalas & Weibel-Mihalas 1984). Figure 16 shows XRoss and Kpianck vs. temperature 
for maximum grain sizes l//m and 1mm. 
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Fig. 16. — Planck and Rosseland mean opacities vs. temperature for a max = 1/mi (solid 
heavy curves) and 1 mm (lighter curves). 



The dips seen in Figure 16 are due to the vaporization temperatures of the different dust 
constituents. For a typical number density of the initial model, about 3 x 10 11 cm -3 , water 
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ice is the main contributor to the opacities under 160 K. At higher temperatures the ice is 
vaporized, so most of the extinction is due to organic grains (< 470 K), troilite (< 740 K), 
and silicates (< 1140 K). Notice that for a max = 1/xm the Planck mean opacities are larger 
than the Rosseland mean opacities in the temperature regime of the simulations (T < 300 
K). This means that for the RC disk, the first term on the right hand side of equation (7) 
is likely to be larger than the second term, so the net cooling of the atmosphere is positive. 
On the other hand, the Rosseland opacities are larger than the Planck opacities for a max = 1 
mm in the same temperature range. In this case, the atmosphere is likely to have negative 
cooling or net heating in equation (7), which, when added to heating by shocks, can cause 
the atmosphere and the disk to expand very quickly and unphysically, as was observed in 
some preliminary tests. To avoid this problem, the calculations presented in Cai et al. (2006) 
use the Planck mean opacities for all terms in equation (7). 

The opacities used here are a function of grain size. Smaller grains contribute more to 
the opacities at temperatures of several hundred Kelvin, while larger grains are the main 
contributors for temperatures lower than 100 K. This indicates that the initial vertical physi- 
cal thickness of the atmosphere in the disk models is determined by the choice of a max . Tests 
show that the atmosphere is physically thinner for a max = 1 mm (higher column optical 
depths) for the temperature ranges of the initial model (tens of Kelvin in the mid and outer 
disk) . This accounts for the longer overall cooling times and slower evolution of the a max = 1 
mm simulation in Cai et al. (2006). The energy transport processes in disk simulations 
depend significantly on the grain size adopted for the opacities. 



As argued here and in Cai et al. (2006), different treatments of boundary conditions 
could drastically change the behavior of the simulation. Although an analytic solution for 
radiative hydrodynamics in a disk does not exist, there are several simple test cases and 
approximations (e.g., Hubeny 1990) that can be explored as a means to test the accuracy 
of any particular radiative physics scheme. We present one such test case in this appendix 
and challenge all researchers who publish radiative hydrodynamics simulations to perform 
similar tests or to develop tests of their own and publish the results. 

The test we present here checks the code's accuracy in achieving a plane-parallel gray 
atmosphere solution. Assume the opacity n is constant, so that 



where m(z) is the surface density down to z. Furthermore, by adopting a constant gravita- 
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tional accleration g, vertical hydrostatic equilibrium requires a pressure profile that follows 

P(z)=gm(z), P(r)=g-. (B2) 

K 

If we require that flux be conserved, then the temperature profile should be close to that of 
the standard Eddington approximation, i.e., 



T(t) 



3 / 



1/4 

(B3) 



where q = 2/3 classically (Eddington approximation) but is closer to q — l/y/S for real 
atmospheres (Mihalas & Weibel-Mihalas 1984). For the test simulation, we place a slab of 
material, which is in hydrostatic, but not radiative, equilibrium vertically in our cylindrical 
grid with a low resolution of (32,8,64). We freeze the radial direction to avoid boundary 
condition problems at the disk edges. In this vein, the test calculation is two-dimensional 
hydrodynamically, but three-dimensional radiatively. The radial flux at the edges of the disk 
is set to zero to conserve vertical flux. A flux crT e 4 ff is introduced in the vertical direction at the 
midplane, where T e s = 800 K for the test. Moreover, we set the ftpianck = XRoss = ac = 1 cm 2 /g 
and m(0) = 20 g/cm 2 so that r m id p iane = 20. 

Figure 17 shows that the temperature structure deviates no more than about 6 % from 
the analytic solution, with q = 2/3, in the region where r > 2/3. The boundary condition we 
impose allows for the correct flux, as calculated by the flux-limited diffusion routine, through 
the disk interior to within a percent despite the error in the temperature. Unfortunately, 
the atmosphere does experience a sudden drop in temperature above the photosphere. This 
drop is due to the lack of the complete cell-to-cell coupling in equation (7) that radiative 
transfer requires. As seen in Figure 17, we expect the RC disk's atmosphere to have a 
temperature drop that is comparable with the drop in this test or even larger due to the 
differences in the Rosseland and Planck means. Because we calculate the boundary flux 
accurately, the temperature drop should not be a problem for disk cooling. The largest 
problem that might be introduced by this sudden drop is dynamical. Recall that the entropy 
S ~ In (P/ p 1 ) ~ In (T/ p 7_1 ), so a sudden temperature drop without a balancing density drop 
will result in a negative vertical entropy gradient. As a result, our atmospheric fitting routine 
apparently seeds superadiabatic regions in our disk where the interior and atmosphere meet 
(Fig. 14). It should be noted that, even so, we do not see thermal convection during the 
Gl-active phase and that we do not expect these superadiabatic regions to result in thermal 
convection because the optical depths are small. 

Our scheme, though crude, accomplishes everything we require: Most importantly, we 
estimate the photospheric flux from the disk interior reasonably well. Although the atmo- 
sphere, which is optically thin, is too cool, we do allow cooling in the atmosphere through 



-41 - 



the first term in equation (7), and we prevent the atmosphere from contracting completely 
onto the disk by adding the second term in equation (7). It is necessary to include the at- 
mosphere in the cooling treatment at some level of approximation because dynamics in this 
region create thermal effects (Boley & Durisen 2006) and because we treat cases with exter- 
nal irradiation (Cai et al. 2006). We are developing additional tests to assess time-dependent 
thermal fluctuations. 




Fig. 17. — Temperature profile for the test simulation. The temperature structure for the 
interior disk (black curve) matches the analytic value (dash curve for q = 2/3) fairly well. 
The atmosphere goes to an isothermal profile, as expected, but it is too cold. The vertical 
line shows the approximate location of r = 2/3. The inset of this figure shows the percent 
change in the flux as calculated by the flux-limited diffusion routine from the Eddington flux 
F e that is introduced at the base of the atmosphere. 
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